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PROGRAMING  LANGUAGE  FOR  THE  SOLUTION  OF  PARTIAL 
DIFFERENTIAL  EQUATIONS  USING  HYBRID  COMPUTERS 

PHASE  I 


l.  INTRODUCTION 

1.  Objective.  The  objective  of  this  report  is  twofold.  The  first  objective  is  to 
provide  a progress  report  on  hybrid-computer  solution  techniques  for  partial  differential 
equations.  The  second  objective  is  to  provide  documented  details  of  the  solution  me- 
chanics and  to  illustrate  the  power  and  speed  of  the  hybrid  computer  when  solving  partial 
differential  equations  (PDE).  A solution -speed  comparison  between  the  hybrid  and  digi- 
tal techniques  shows  the  hybrid  to  be  the  faster  of  the  two. 

2.  Background.  The  Electrical  Equipment  Division,  U.S.  Army  Mobility  Equip- 
ment Research  and  Development  Center  (USAMERDC),  is  involved  in  the  research, 
development,  and  engineering  of  electromagnetic  machinery,  power  conditioners,  and 
power  electronics  components  (SCR’s,  transistors,  and  rectifiers).  These  efforts  require 
the  solution  of  partial  differential  equations  in  order  to  provide  flux  plots  and  equipo- 
tential  plots.  When  digital-computer  techniques  arc  used,  these  problem  solutions  arc 
slow  and  costly.  However,  by  using  hybrid-computer  techniques,  we  can  reduce  these 
computing  costs  by  a factor  of  15  to  25,  with  a corresponding  increase  in  computing 
speed  by  a factor  of  between  15  and  100.  The  Electrical  Equipment  Division  has  a 
powerful,  interactive  hybrid-computer  facility  (Figure  1),  which  is  part  of  the  CAD-E 
facility  (Figure  2).  The  hybrid  computer  is  a Digital  Equipment  Corporation  PDP-15/ 
Applied  Dynamics  AD-4  hybrid  computer  coupled  to  a Tektronix  4010  Graphic  Termi- 
nal. Figure  3 shows  the  PDP-15/76  digital  processor  which  has  a unichannel,  1.2- 
million-word  disk  and  16K  of  core.  The  AD-4  analog  processor  (Figure  4)  has  96  ampli- 
fiers as  well  as  an  autopatch  capability.  The  technical  paper  Hybrid  Computer  Solution 
Techniques  for  Laplace’s  Equations,  by  the  authors  of  this  report,  has  helped  immensely 
in  preparing  this  report.* 

3.  Organization.  This  report  is  divided  into  five  parts:  Introduction,  Program 
Philosophy,  Computer-Solution  Mechanics,  Examples,  and  Conclusions  and  Future  Work. 
Additional  material  is  given  in  the  three  appendixes.  The  Program  Philosophy  section 
describes  the  philosophy  of  program  development.  The  section  on  Computer-Solution 
Mechanics  presents  the  details  of  problem  setup  for  the  hybrid-computer  solution.  The 
Examples  section  and  the  appendixes  present  sample  problem  solutions  and  special  con- 
siderations. This  report  will  provide  the  basis  for  comparing  the  interactive  hybrid- 
computer  solution  of  partial  differential  equations  to  the  digital-computer  approach. 

* J.  T.  Broach  and  H.  M.  McKechnic,  Hybrid  Computer  Solution  Techniques  for  Laplace’s  Equation,  Proceedings  of 
1974  Anny  Numerical  Analysis  Conference,  AUO  Report  74-2,  pp.  253-271. 


1 


PDP15  DIGITAL  PROCESSOR 


AD-4  ANALOG  PROCESSOR 


TEKTRONIX  4010  GRAPHICS  TERMINAL 


Figure  1.  Interactive  hybrid-computer  facility. 
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Figure  2.  Computer-aided  design  engineering  facility. 
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II.  PROGRAM  PHILOSOPHY 


This  report  describes  a hybrid-computer  solution  approach  to  the  solution  of 
partial  differential  equations.  However,  to  understand  the  reasoning  for  this  method, 
the  pure  analog-computer  approach  to  the  solution  of  partial  differential  equations 
must  he  discussed.  The  technical  background  for  this  effort  also  will  be  useful  to  the 
full  understanding  of  the  program  philosophy. 

4.  Technical  Background.  The  background  of  the  present  work,  typical  equa- 
tions, and  their  method  of  solution  are  discussed  below. 

a.  Status  in  this  Area  of  Work.  During  the  early  1960’s,  much  work  was 
accomplished  for  the  solution  of  partial  differential  equations  on  analog  computers. 

With  the  expected  use  of  hybrid  computers,  the  emphasis  was  shifted  to  their  utilization. 
However,  the  efforts  since  then  have  been  small,  with  little  to  show  but  theory.  In  the 
digital  area,  work  has  progressed,  mainly  because  of  the  easier  man/machine  interface 
and  because  of  the  efforts  of  universities  and  the  large  computer  companies. 

b.  Types  of  Problems.  The  Electrical  Equipment  Division  is  involved  in  the 
solution  of  partial  differential  equations  for  heat  transfer  and  magnetic  flux  in  electric 
and  electronic  equipment.  As  a result,  the  first  problem  to  be  examined  and  set  up  will 
be  the  diffusion  problem  and  its  associated  equations.  The  solution  of  this  type  of  equa- 
tion will  provide  immediate  benefits  to  the  Electrical  Equipment  Division. 

c.  Types  of  Partial  Differential  Equations.  There  are  three  types  of  partial 
differential  equations  which  are  representative  of  a large  number  of  engineering  problems 
encountered: 


K — = V2  0 + f 
at 

(heat  equation  or  diffusion  equation), 

=v2  0 + f 

at2 

(wave  equation),  and 

K^  = V%  + f 
at2 

(dynamic  structural  equation  (biharmonic 
equation)), 

whcreV%  = + |!|  + |!i 

9x2  9y2  9z2 

and V4  ^ = ?A  + ^+*l±  + 
0x4  3y4  flz4 

+ „ I a4*  „ + a4*  l 

| a*2  3y2  ax2  a*2  ay2  a*2  f 
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d.  Usual  Methods  of  Solution.  There  are  three  major  techniques  of  solu- 
tion: (1)  separation  of  variables,  (2)  finite  difference,  and  (3)  stochastic.  Generally, 
we  will  use  the  finite-difference  technique  because  it  can  handle  time-varying  boundary 
conditions  and  nonlinearities  easily.  The  separation-of-variables  technique  assumes 
linearity.  For  the  digital  solution,  one  reduees  the  partial  differential  equation  to  a set 
of  algebraic  equations  using  the  finite-difference  technique.  This  means  that  iterative 
techniques  must  be  employed  to  obtain  solutions.  For  the  analog  solution,  one  obtains 
a set  of  ordinary  differential  equations  using  the  finite-difference  technique. 


5.  Analog  Approach.  The  general  approach  to  be  used  to  solve  the  two- 

dimensional  Laplace  equation,  = 0?  is  to  use  finite  differences  for  one  of 

the  space  variables  and  to  solve  the  other  variable  continuously.  On  the  analog  compu- 
ter, this  means  we  have  two  choices.  We  can  divide  the  space  such  that  we  solve  for  a 
continuous  solution  as  a function  of  y at  each  of  a series  of  x-stations,  or  we  can  solve 
for  a continuous  solution  as  a function  of  x at  each  of  a series  of  y-stations.  Basing 
our  calculations  on  engineering  considerations  for  accuracy,  we  will  try  to  use  only  a 
few  stations.  This  also  will  reduce  the  number  of  analog  components.  In  order  to 
demonstrate  solution  accuracy  and  to  identify  mechanization  problems,  the  first  test 
problem  is  one  that  has  an  exact  solution  and  that  is  a special  ease  of  the  more  general 
problem  which  will  be  solved  as  the  approach  is  refined  into  a programing  language. 


The  interesting  general  problem  for  the  electrical  engineer  designing  military 
generators  and  motors  is  one  which  provides  the  flux  or  flowline  patterns  and  the  equi- 
potential-line  patterns  of  the  magnetostatic  field  in  a section  of  the  air  gap  of  the 
machine.  Figure  5 is  a diagram  of  this  complicated  geometry.  Here  we  need  to  be  able 
to  take  care  of  a complicated  geometry  with  different  types  of  iron  and  with  various 
boundary  conditions.  The  overall  objective  is  to  provide  a language  which  allows  the 
design  engineer  to  draw  this  picture  on  the  graphic  screen,  to  input  the  required  bound- 
ary conditions,  to  solve  the  problem  on  the  hybrid  computer,  and  to  provide  a picture 
of  the  desired  distributions  of  flux  and  potential,  displayed  on  the  graphic  screen.  The 
first  test  ease  is  a simplified  example,  that  will  allow  for  an  exact  solution,  which  can 
be  used  for  a comparison  of  results.  Figure  6 is  a diagram  of  a rectangular  space  used 
for  the  first  test  case. 


III.  COMPUTER-SOLUTION  MECHANICS 


6.  Solution  Mechanics.  For  the  test  ease,  we  have  a rectangular  region,  and  we 
will  investigate  the  field  inside  this  region  when  three  boundaries  arc  at  \p-0  and  one  is 

at  \//-f(x).  The  exact  solution  for  this  case  is  100  \[/(x,y)=100  sin  — . ^ - , 

a sinh  [ 7r  b/a  ] 
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A B 


Figure  5.  Typical  electromagnetic  machine  geometry. 
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where  a and  b arc  as  defined  in  Figure  6.  This  solution  has  been  mechanized  on  the 
PDP-15  section  of  the  hybrid  to  provide  \f/(x,y)  for  comparisons.  Two  analog  solutions 
have  been  studied  (6a  and  6b,  below). 


a.  Continuous  x,  Discrete  y.  In  this  solution,  the  analog  computer  simul- 
taneously solves  a set  of  differential  equations  at  each  of  a series  of  y-stations  to  pro- 
vide ip(x)  |y  , where  a is  the  station  number/location,  which  will  give  the  value  of  i//(x, y) 
at  all  points  if  y is  on  a station  line.  Some  extrapolation  means  is  assumed:  of  course, 
if  Ay  is  small  enough,  it  will  not  matter.  For  this  solution-method  example,  we  will  use 
six  stations  in  the  y-direction.  For  boundary  conditions,  we  have  even  derivates  (y0  is 

7TX 

considered  even)  specified  at  the  boundaries  y = 100  sin  — , and  y5=  0 for  all  x.  Also, 

° a 

y,,y2,y3,  and  y4  have  a boundary  condition  of  0 for  x=0  and  x=a.  In  Hausner’s rules 
for  mechanization  (Appendix  A),  rule  2 states  that  we  should  arrange  the  grid  station 
so  that  an  integer  station  (yo,  ys)  appears  at  the  boundary  since  we  have  even  derivatives 
specified  at  the  boundary.  The  next  Hausner  rule  (rule  3)  says  that  we  should  generate 
high-order  derivates  with  first-order  approximations,  mechanizing  all  lower  order  de- 
rivates as  summational  outputs. 


Thus,  we  let  IT  = \I'” 


*j.l  - 2 + 


where  h is  and  j is  » 


h2 

, so  Dj  « 


and  ^ 


- \b. , + i b. 
e.  i-1  LI 


- V + 


J+>4 


Thus,  we  gene- 


h J h 

rate  five  intermediate  solutions  (0j/2’  ^3/2’  ^5/2’  ^7/2’  8,13  09/2)  an(^  use  e*ght  integra- 
tors (Figure  7). 


Setting 


31 


9y2 


32^n 

3x2 


3y2 


equation  gives  us: 


n +Yi 

(Ay)= 

31 'fc 


3x2 


, a finite-difference  equation  for  y,  in  the 


Then  we  can  solve 


,_l 

yn 

(Ay)2  J 

for  \p  (x)  |y  by  using  the  unsealed  equations: 


d2\!/, 

_ ^1/2“  ^3/2 

dx2 

(Ay)2 

d2^ 

_ ^3/2“  ^5/2 

dx2 

(Ay)2 

<l20s 

_ ^5/2"  4*7/2 

dx2 

(Ay)2 
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mY)--o 


TfCL.Y)  =0 


^'igu^c  6.  Rectangular  space. 
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Figure  7.  Grid  for  continuous  x,  discrete  y. 
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d2<//4  _ 07/2-<ft9/2 


dx2 

(Ay)2 

^1/2  = 

^1 

- 100  sin  — 
a 

^3/2  = 

'Pi 

- 'Pi 

^5/2  = 

\ l>  3 

-'Pi 

^7/2  = 

'Pa 

- 'Pi 

&)/2  = 

0- 

'Pa 

For  mechanization  purposes,  we  replace  t by  x in  a one-to-one  replacement  (i.e.,  1 
second  = 1 unit  of  distance  in  x). 


In  the  unsealed  equation,  Ay  = — , so  we  have  a way 

4 y (No.  of  Stations -1)  y 

to  incorporate  a and  b in  the  solution.  For  scaling  use  the  values  given  in  the  table  are 

typical. 


Variable 

Est.  Max.  Value 

Scale  Factor 

Scaled  Computer  Variable 

0 

100  v 

100 

100 

'P 

100  v 

100 

100 

m 

'p1 

100  v/s 

100 

100 

m 

(In  our  problem  as  it  is  set  up,  the  X-generator  (Integrator  271)  is  generating  10  v/s  or 
0.1  s/v.  When  we  measure  10  volts  on  X at  10  v/s,  we  had  1 second,  or  1 unit  of  dis- 
tance in  X,  which  corresponds  to  a.)  For  tliis  problem,  we  used  the  initial-condition 
(IC)  pots  on  the  ^'-integrator  to  obtain  the  proper  boundary  condition  for  through 
\pA  at  x=a.  In  this  problem,  we  used  these  pots  to  make  \[/i , ip2,  \p3,  and  ipA=0  at  x=a. 

b.  Continuous  y,  Discrete  x.  This  method  is  identical  to  the  continuous  x, 
discrete  y method  except  that  the  problem  space  is  divided  into  stations  in  the  x-direction. 
The  problem  is  solved  continuously  in  the  y-dircction.  This  method  is  discussed  in  more 
detail  in  the  examples  (section  IV). 
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7.  Special  Techniques.  Two  special  techniques  for  problem  solution  may  be 
mentioned. 


a.  Dividing  Problem  Space.  In  an  effort  to  minimize  equipment  and  to 
provide  an  easy  conversion  to  autopatch,  we  will  divide  the  problem  space  into  three 
fixed  stations  and  one  variable  station.  Using  symmetry  (special  case),  we  get  mirror- 
image  solutions  in  the  right  half  and  in  the  left  half  of  the  rectangular  space.  Therefore, 
by  this  consideration,  we  get  2n-3  solutions  for  n stations.  Using  the  hybrid-solution 
control,  we  will  set  the  variable  station  at  a specified  AX -spacing  from  the  center  station, 
and  a solution  will  be  obtained.  Then  AX  will  be  increased,  and  the  problem  will  be 
solved  again.  This  iterative  process  will  be  repeated  until  all  specified  stations  are  used. 
This  method  allows  for  linear  or  nonlinear  spacing. 

b.  Approaching  Boundary  Value  by  Varying  IC-Pots.  Another  iterative 
process  found  to  be  useful  occurs  in  satisfying  the  boundary  equations.  By  varying  the 
IC-pots  on  the  ^/-integrators  one  at  a time  and  in  station  order  from  left  to  right,  we  can 
iteratively  approach  the  required  boundary  value.  This  method  requires  that  the  first 
pot  1m*.  varied  until  the  \J/j  -variable  equals  zero  at  the  prescribed  location  on  the  x-axis 
(x=a)  while  all  other  pots  are  fixed.  Then,  the  second  pot  is  varied  until  \p= 0 at  the 
same  location.  This  process  is  repeated  sequentially  until  all  variables  (rpi , i//2,  ^3, 
and  \]/A)  are  zero  at  the  same  point.  This  method  will  be  illustrated  clearly  by  the 
examples,  which  follow  in  the  next  section.  Both  of  the  iterative  processes  described 
above  arc  performed  rapidly  by  the  PDP-15  digital  computer. 

IV.  EXAMPLES 


8.  Laplace  Equations  for  Two-Dimensional  Solution.  The  geometry  of  this  prob- 
lem dictates  use  of  the  continuous  y,  discrete  x solution  method.  Based  on  trial  solu- 
tions, it  was  determined  that  six  stations  are  adequate  (five  fixed  and  one  variable  station). 
Two  stations  are  at  the  boundaries,  x=0  and  x=a,  where  \f/0=  i//5=0.  Figure  8 is  a dia- 
gram of  the  space,  with  the  variable  station  shown  as  a broken  line. 


For  this  mechanization,  XQ,  Xi , X2,  X3,  and  X5  are  fixed  locations,  and 
X4  varies.  Because  of  symmetry,  Xj  and  X2  will  have  mirror-image  solutions  in  the 
right  half-space,  and  X*  will  have  mirror-image  solutions  in  the  left  half-space.  Point 

2sl  • • 

X3  is  located  at  — , while  X,  is  at  a/6  and  X2  is  at  — . By  symmetry  conditions, 

6 6 

there  will  be  an  identical  solution  to  X2  at  — , to  Xi  at  ~ , and  to  X4  at  — + — , 

6 6 6 


with  K4  being  specified  by  the  user.  For  initial  conditions  along  the  y=0  boundary. 
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X» 


X 


Figure  8.  Problem  space. 


!//o=  100sin(^  (0))  ; = 100  sin  (1  (a/6))  ;^=100sin  0 ; ^3=100 

Sin  (a  (^))  ;*»  = 100sin  (f  (f ))? *4=  100sin(l  (^))  ;Where 
K = 3 + K4 . 

When  the  method  described  previously  was  used,  it  was  possible  to  solve  the 

equations: 


Equation 

Definition 

= 

(^l/2~  ^3/2) 

Ax,, 

a a 

Ax„  -- 

(i) 

4>'i 

= 

(03/2“  ^5/2) 
Ax2| 

Ax2,  =^- 
6 

(2) 

*3 

(05/2r  ^7/2) 
Ax3, 

AX3,  ='^(|)  + V4(K4) 

(3) 

*4 

= 

(^7/2“  ^9/2) 
Ax4, 

Ax4,  ='/4(K.)+S4 

(i-K‘) 

(4) 

^1/2 

= 

(vt/i-  \fO 
Ax, 2 

Ax,^| 

(5) 

^3/2 

= 

A Xj2 

AxI2  =| 

(6) 

^5/2 

= 

('I'i-'l'l) 

Ax32 

a a 

Ax32-- 

(7) 

^7/2 

= 

(^4-^3) 

Ax42 

A X42  = K4 

(8) 

^9/2 

= 

(i/'s-'M 

Ax52 

A x52  = — - K4 
0 

(9) 

Variable  K4  is  defined  as  follows: 

K4=Kr(|2)  , (10) 


where  KR  is  the  spacing  factor. 
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Changing  the  equation  form,  we  obtain  the  following  \p  and  <j>  values: 


$1  = ( 

(*l/2_  *3/2  ^ 

(11) 

$2  = ( 

— ) 
^ A X^j/ 

(*3/2-  ^5/2) 

(12) 

4>3  = ( 

<*5/2~  *7/2> 

(13) 

4>4  = ( 

a 

<*7/2  — *9/2) 

(14) 

*J/2  = ( 

^Ax,2) 

Wl“  *0) 

(15) 

*3/2  = ( 

' 1 } 
^ AXj2  / 

Wa-  ^1) 

(16) 

*5/2  = ( 

it) 

(^3-  '/'a) 

(17) 

*7/2  = ( 

1 ) 
Ax42  / 

(^4-  '/'a) 

(18) 

*9/2  = ( 

1 ) 
Ax52  / 

(^s-  <M 

(19) 

Continuing  to  change  the  equation  form  (since  \po=- 

following: 

\[/5=0),  we  obtain  the 

(A  x,2)</> 

1/2  ~ *» 

(20) 

(A  x22)0 

3/2  = *2 

- ^1 

(21) 

(Ax32)  </> 

5/2  = *3 

- V'a 

(22) 

(A  x42 ) <f> 

7/2  = *4 

- ^3 

(23) 

(A  xS2)</> 

9/2  = <" 

(24) 
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001  +■  = (Hr,)  ♦.»  • (Hr,)  *w  <25) 

°-01*’  = (Hr,)  ^-(Hr,W  (26) 

001  ■ (Hr,)  ♦w-(i^)  *«•  <27) 

001  *•"  - (Hr,)  ♦»*  - (Hr)  *«  <28> 

(Axij)0i/2  = (P224)  0i  (29) 

(AX22)  03/2  = (P226)  02  - (P223)  0i  (30) 

(Ax32)05^2  = (P244)  03“  (P236)  02  (31) 

(K,Ax42)07/2  =(K1)(Pj66)  04-(K1)(P247)  03  (32) 

K,  = ^1  (33) 

Ax42 

(K2  AxS2)09/2  = -(K2  P256)  04  (34) 

K2=^21  (35) 

AxS2 

Continuing  the  rearrangements: 

1*232  (Ax,2)0]/2  = (j“)  01/2  (36) 

P245  (AX22)  03/2  = (^~)  03/2  ^37) 

P227  (Ax22)03/2  = 03/2 

P233  (Ax32  ) 05/2  = (l^i)  05/2  (39) 

P243  (AX32)  05/2  = (2^7.)  ^5/2 
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Ki  P253  (^X42)^7/2  = 

0 x1 

ol< 

^7/2 

(41) 

• 

Ki  P265  (^x42)^7/2  = 

\A  x4|  / 

^7/2 

(42) 

^2  P276  (^x52  ) ^9/2  = 

(™l) 

\Ax41/ 

<A9/2 

(43) 

Finally,  we  obtain  the  following  pot  settings: 

P224  = 1 

(44) 

P226  = 1 

(45) 

P223  = 1 

(46) 

P244  = 1 

(47) 

P236  = 1 

(48) 

P _Ax32 

* 266  A” 
Ax42 

(49) 

• 

p -Ax32 

r247  A 

Ax42 

(50) 

p — ^x32 

1 256  A 

AxS2 

(51) 

p _ 0.01 

(52) 

232  (Axn)(Ax12) 

p _ 0.01 

(53) 

(Axn)  (Ax22) 

p _ 0.01 

(54) 

(A X22 ) (Ax2j  ) 

p _ 0.01 

(55) 

(Ax32)  (Ax2|  ) 

p _ 0.01 

(56) 

243  (Ax32)(Ax31) 
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• 

(57) 


(Ax3! ) (Ax42)  (Kj  ) 

0.01 

(Ax4i)  (Ax42)  (Ki) 


(58) 


1*276 


aoi 

(Ax41)  (Axs2)  (k2) 


(59) 


The  program  will  scan  the  space  as  previously  described,  and  with  four  differ- 
ent positions  for  station  X4  we  actually  obtain  data  for  15  equivalent  stations  as  is  shown 
by  Figure  9. 

In  order  to  obtain  the  desired  plots,  it  is  necessary  to  perform  a core  search 
for  a specified  (//-value: 

a.  Check  out  the  specified  X-station  and  its  equivalent  image. 

b.  Use  straight-line  interpolation  between  data  points. 

For  example:  y value  = ITM/10,000,  where  ITM  = b 

x-value  = x-station  location 


For  a specified  X: 

a.  Start  at  the  maximum  ^-value  until  ^ in  core  is  less  than  the  specified 

b.  Back  up  one  space  and  check  discrete  y-values;  use  linear  extrapolation 
to  get  specified  value  x,y  data. 

9.  Hybrid-Computer  Solution. 

a.  General.  The  hybrid-computer  solution  may  be  illustrated  graphically. 
The  problem-space  geometry  is  shown  in  Figure  8 and  the  space  with  the  solution  grid 
is  shown  in  Figure  9.  The  finite-difference  equations  are  shown  in  Figure  10,  and  the 
computer  patching  diagram  is  given  as  Figure  11.  A program  control  flow  chart  is  shown 
in  Figure  12,  and  the  patchboard  is  shown  in  Figure  13.  Figure  14  shows  the  logic  patch- 
board. 


b.  Computer  Program  PDR2B.  The  computer  program  is  stored  in  the 
execute  file,  PDR2B,  in  the  RMM  file  on  disk.  Program  listings  and  subroutines  are 
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BASIC  FINITE  DIFFERENCE  SCHEME 
FOR  HYBRID  COMPUTER 


A CHANGE  TO  AN  ORDINARY  2nd  ORDER  DIFFERENTIAL 
EQUATION  AT  EACH  X-STATION 


♦ = 


1 AX 


11 


i+w- 


* 1 WHERE 

' 3/2  1 


to 


♦ = 


2 AX 


I * 3/2  *s/9  I 
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t = 


3 AX 


31 


1^5/2"  * 


5/2 


7/2 


AX 


) I* 


41 


7/2 


- d>  I 

'9/2  1 


4>  = | — ■ — I I \i/  — \i/ 

M/2  lAX12MY1  Y 


1 


*3/2  = [Zlf721  1 +2~  + l 

l 


4> 


5/2  'AX 


1 


I I*,"* 


32 


4v,  = (tt-I  I'K-'J' 


7/2  'AX 


42 


*9/2  ^AX5J  ^5  *4 


Figure  10.  Finite-difference  equations. 


L-OO 


■oo — 


Figure  11.  Computer  patching  diagram. 


INPUT: 


USER 


3 


RECTANGLE  BOUNDARIES 
# OF  GRID-LINES 


TEK  4010 


PDP-15  DIGITAL  PROCESSOR 


AD-4  ANALOG  PROCESSOR 


PDP-15  DIGITAL  PROCESSOR 


TEK  4010 


Figure  12.  Program  flow  chart. 
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Figure  13.  Analog  patchboard. 
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Figure  14.  Logic  patchboard. 


given  in  Appendix  B.  The  large  size  of  this  problem  requires  “chaining,”  and  the  pro- 
gram details  are  in  Appendix  C.  The  following  is  a description  of  the  use  and  response 
of  PDR2B.  With  the  PDP-15/AD-4  hybrid  up  and  running,  the  PDP-15  executive 
supplies  a “$”  to  indicate  user  input.  To  the  “$”  on  the  Tek»'  mix  4010,  the  user 
types  in  “E  PDR2B.”  The  computer  prompting  response  is  a statement  for  input: 
“Input  A,  B,  DEL1,DEL2,  DEL3,  DEL4,  DEL5:  F5.2,  5F5.4.”  This  allows  the  user 
to  provide  the  x and  y space  dimensions  (A  and  B,  respectively).  The  spacing  for  the 
variable  grid  line,  referenced  from  the  center  line,  is  not  used.  Once  this  spacing  is  in- 
put, the  computer  responds  with  the  prompting:  “Specify  Number  of  Lines  LT16.” 
Tills  allows  the  user  to  vary  the  number  of  stations  for  trial  solutions.  The  computer 
prints  the  value  of  DEL  (as  measured  from  the  center  x-station)  and  the  IC-pot  values, 
which  are  required  to  satisfy  the  boundary  conditions  through  the  closed-loop,  analog 
iterative  process,  described  in  Appendix  B.  Figure  15  shows  ihe  computer  prompting. 
Program  solution  output  is  shown  by  Figures  16  and  17.  Figure  18  illustrates  the  solu- 
tion with  a grid,  while  Figure  19  depicts  the  solution  without  a grid.  Normally,  for  pro- 
duction runs,  the  problem  grid  would  be  well  specified;  but,  for  this  problem,  it  was 
not.  Several  linear  and  nonlinear  spacings  were  investigated.  It  should  be  noted  that 
the  nonlinear  grid  helps  to  clarify  solution  slopes  in  specific  areas  of  interest.  The  use 
of  nonlinear  grid  is  optional  (i.c.,  it  can  be  selected  as  nec.ded).  The  16K  core  of  the 
present  PDP-15  digital  subsection  of  the  hybrid  unit  limits  us  to  about  20  grid  stations 
(40  with  symmetry),  but  more  would  be  available  if  we  had  written  the  solution  to 
disk  or  tape  storage  and  had  performed  the  graphics  with  another  program.  Also,  the 
graphics  display  uses  a simplified,  point-to-point  plotting  routine,  which  could  be  re- 
fined for  smoother  curves. 

The  b/a-ratio  limits  for  this  method  as  it  is  presently  programed  are 
between  0.1  and  0.3,  mainly  because  of  the  assumed  scaling.  This  limitation  will  be 
eliminated  later,  but  it  is  not  serious  enough  to  warrant  a change  for  the  trial  example. 
Figures  15  through  19,  which  depict  the  solution  on  the  Tektronix  4010  Graphic  Ter- 
minal screen,  were  used  to  demonstrate  the  problem  I/O  and  do  not  describe  accurate 
solutions.  The  next  set  of  figures,  which  is  hardcopy  output  for  the  Tektronix  graphics 
display,  is  used  to  provide  the  comparison  of  accuracy  between  the  exact  and  hybrid 
solutions  for  this  example.  The  exact  solution  uses  a mathematical  solution  subroutine 
in  place  of  the  hybrid  subroutine  set  PDE,  MCON,  and  PDE2  (see  Appendix  C for  more 
details).  All  other  input  and  output  subroutines  stay  the  same.  Using  the  problem  def- 
inition parameters  (A=l,  B=.l)  and  10  lines  (stations),  wc  can  compare  results.  Note 
that  the  computer  uses  nine  lines  to  divide  the  right-hand  space  of  the  problem  into  10 
spaces.  Figure  20  is  the  hardcopy  output  for  the  hybrid  solution,  and  Figure  21  is  the 
hardcopy  output  for  the  exact  solution.  Appendix  C contains  \!/(y)-data  for  each 
X-station  generated  by  the  exact  and  hybrid  solutions. 
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Figure  17.  Program  solution  output  (completed). 
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Figure  18.  Solution  with  grid. 
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Figure  20.  Hardcopy  output  of  hybrid  solution. 


Figure  21.  Hardcopy  output  of  exact  solution. 
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In  clock  time,  each  hybrid-computer  solution  set  took  30  seconds.  (A 
33-grid  solution,  including  the  symmetry,  took  about  7 minutes.)  The  hybrid-computer 
solution  runs  100  times  faster  than  real  time  and  is  faster  than  the  exact  solution  pro- 
vided by  the  PDP-15  only.  Figures  16  and  17  verify  our  original  assumption:  that  we 
could  scan  the  space,  while  maintaining  five  stations  fixed  and  one  moving,  because 
the  first  three  pot  settings  (two  stations  are  at  the  boundary,  where  \b=0)  always  return 
to  the  same  value  at  solution;  however,  the  grid  station,  being  moved,  changes  the  pot 
value. 

V.  CONCLUSIONS  AND  FUTURE  WORK 

10.  Conclusions  and  Future  Work.  So  far,  we  have  shown  a technique  for  solving 
partial  differential  equations  on  hybrid  computers  which  is  at  least  50  times  faster  than 
the  digital  solution.  This  speed  of  solution  occurs  because  we  solve  the  problem  in  a 
continuous,  closed-loop,  analog  process.  Also,  we  have  established  an  iterative  solution 
technique,  which  converges  rapidly  and  allows  us  to  maintain  overall,  simplified  digital 
control  over  the  closed-loop,  analog  solution  process.  The  comparison  of  the  hybrid 
solution  to  the  exact  analytical  solution  demonstrates  the  accuracy  of  this  approach. 

The  next  steps  are  to  generate  the  problem  menus  and  to  solve  the  field 
problem  for  a slot  geometry  and,  then,  for  other  complex  geometries.  The  progress 
demonstrated  to  date  offers  an  optimistic  outlook  for  complete  success  in  the  future 
planned  work  of  this  project. 
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APPENDIX  A 


HAUSNER’S*  RULES  FOR  MECHANIZATION 


The  following  is  a list  of  Hausner’s  Rules  used  in  this  project: 

Rule  1 — To  obtain  a kth-ordcr  solution,  all  approximations  must  be  kth 
order,  including  those  accounting  for  boundary  conditions. 

Rule  2 — If  only  even  derivatives  of  a dependent  variable  (such  as  u,  u",  u"  ", 
etc.)  are  specified  at  a boundary,  arrange  the  grid  stations  so  that  an  integer  station 
(say,  XD  or  Xj ) appears  at  a boundary.  If  at  least  one  odd  derivative  (u',  u"\  etc.)  is 
specified  at  a boundary,  a half-integer  station  (say,  X^2)  should  be  placed  at  a boundary. 

Rule  3 — Generate  high-order  dcrivates  with  first-order-derivative  approxima- 
tions, mechanizing  all  lower  order  derivatives  as  summational  outputs. 


* 

A.  Hausner,  “Analog  and  Analog/Hybrid  Computer  Programing,”  Prcntice-Hall  1971,  pp  435-436. 
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APPENDIX  B 


ANALOG  CONTROL  ROUTINES 


A brief  discussion  of  the  analog  control  routines  used  to  reach  solutions  is  given 
in  this  appendix. 

B-l.  Differentiation  with  Respect  to  y.  The  analog  computer  actually  performs 

d^  aS  It  ’ w*iere  y ‘s  rcPrescnted  as  t on  a one-to-one  basis.  The  time-base  (or  y-base) 

generator,  integrator  271,  normally  is  providing  10  v/s;  thus,  wc  get  0.1  s/v  as  the  out- 
put. Since  1 unit  of  y is  equivalent  to  1 second,  it  takes  0.2  second  to  provide  0.2 
unit  of  y.  This  means  that  the  integrator  output  is  2 volts  in  0.2  second  (0.1  s/v  * 2 
volts  = 0.2  second).  In  order  to  provide  the  proper  output  rate  for  integrator  271,  pot 
273  is  set  to  0.01  with  100  volts  input.  The  normal  integrator  rate  is  10  v/s  in  quad- 
rant two  of  the  analog  patchboard. 

B-2.  Closed-Loop  Analog  Solution.  The  fastest  possible  solution  is  obtained 
when  the  analog  computer  operates  in  a closed-loop  fashion.  The  solution  control  is 
accomplished  as  follows:  (1)  The  user  provides  input  parameters  to  the  digital  unit; 

(2)  the  digital  computer  uses  these  parameters  to  automatically  scale  the  problem,  to 
set  the  analog  comparator  pot  settings  for  time  (or  b)  value  in  order  to  place  the  com- 
puter in  hold,  and  to  set  the  pots  and  start  the  solution;  (3)  the  digital  unit  waits  a 
sufficient  time  in  order  to  allow  the  analog  unit  to  go  to  “hold,”  checks  end-point 
values  for  convergence,  resets  the  computer  to  run  again,  and  repeats  this  until  con- 
vergence occurs;  (4)  once  convergence  occurs,  the  digital  unit  resets  the  computer  and 
causes  the  analog  unit  to  operate  for  a set  number  of  predetermined  increments,  at 
which  points  the  analog  comparator  places  the  computer  into  the  “hold”  mode  and  the 
digital  unit  samples  and  stores  v//-,  y-,  and  x-data;  (5)  this  process  is  repeated  until  all 
specified  x-stations  have  been  used;  (6)  once  all  x-stations  have  been  used,  the  digital 
unit  asks  the  user  to  specify  \//,  A and  the  number  of  lines  to  plot;  and  (7)  the  digi- 
tal unit  uses  these  data  to  search  its  stored  i //-,  y-,  and  x-data  and  to  provide  the  plot. 
The  digital  unit  is  programed  to  provide  many  variations  of  the  plotting,  once  the  hy- 
brid unit  has  finished  computing,  in  order  to  keep  from  having  to  recompute  each  time 
a new  plot  variation  is  needed. 

B-3.  Analog  Comparator  Logic.  The  logic  and  analog  patching  needed  to  accom- 
plish the  time  (or  y-)  control  is  shown  by  Figure  11.  The  output  of  integrator  271  is 
fed  through  pot  277  to  amplifier  233.  The  output  of  amplifier  233  goes  to  comparator 
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231  on  the  analog  patchboard.  The  reference  voltage  (equivalent  to  y=b)  comes  from 
amplifier  223,  which  is  the  other  input  to  comparator  231.  When  the  sum  of  the  in- 
puts goes  positive  (occurs  at  the  instant  y becomes  infinitesimally  larger  than  b),  a logic 
1 is  generated  by  the  out-point  on  the  logic  patchboard.  Since  “out”  on  comparator 
231  is  connected  to  SYS  Hold,  it  receives  a logic  1,  which  places  the  analog  unit  in  the 
“hold”  mode,  thus  stopping  computation.  In  order  to  reset  properly,  the  digital  unit 
overrides  the  patched  “hold”  mode  by  a “hold”  command,  reads  the  desired  Rvalue, 
places  the  computer  in  the  IC-mode,  and  resets  the  comparator  output  to  logic  0 by 
setting  pot  237  to  0.  For  the  sampling  of  i //-,  x-,  and  y-data  after  convergence  tests  are 
met,  pot  237  is  incremented  to  the  preset  values,  thus  stopping  the  computer  at  the 
desired  points  of  y,  reading  the  data,  and  continuing  to  the  next  point  as  soon  as  pot 
237  is  updated.  This  process  is  limited  to  12  data  points  because  of  the  dimension 
statement,  which  reflects  present  core  limits.  Methods  that  would  allow  more  points 
could  he  used  but  are  not  required  for  the  test  example. 

B-4.  Iteration  Control.  The  computer  is  programed  to  set  all  four  IC-pots  (for 
the  first  derivative  of  i p)  initially  and,  then,  to  go  tlirough  a preset  sequence  to  set  the 
first  IC-pot  until  goes  to  0 at  y=b.  The  computer  then  goes  to  the  second  IC-pot 
and  changes  it  until  \!/2=0  at  y=l>.  Each  time,  all  \Jj'a  are  sampled  to  see  if  they  are 
simultaneously  0 at  y=b.  This  process  continues  to  ^3,  \I/4,  , i^2,  i//3, 1^4,  etc.  until 

\Jjl={[/j=iJ,3=ip4=0  at  y=b.  This  process  generally  converges  in  less  than  30  seconds 
(about  10  iterations  at  most). 
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APPENDIX  C 


COMPUTER  PROGRAM  LISTINGS 


C-l.  Introduction.  This  appendix  gives  a listing  of  the  problem  source  code,  the 
chaining  routine,  and  the  several  programs  used  in  this  study  and  depicts  the  program 
flowcharts. 

The  hybrid-computer  program  consists  of  a main  program  (designated  sub- 
routine POT)  and  eight  subroutines:  PDE,  MCON,  CON,  READSI,  PDE2,  DISK,  DRW, 
and  DRWA.  The  hybrid-computer  program  also  requires  the  hybrid  routines  and  the 
Tektronix  routines.  The  problem  requires  “chaining”  on  the  16K  core  configuration 
of  the  PDP-15.  The  chaining  routine  produces  the  XCT  and  XCU  files  and  allows  the 
program  to  be  run  by  using  E PDR2B. 

C-2.  Hybrid  Program  Listings.  The  following  listings  are  the  routines  used  for  the 
hybrid-computer  solution. 
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c 

c 


2051 

6604 


2650 


601 

900 

69? 


11 


SUBROUTINE  POT  ' X 

THIS  WILL  ACT  AS  THE  MAIN  PROGRAM 
DIMENSION  TST< 2) 

C0MM0N/UI/YC18),  I PS  I 

COMMON/BRWP/XLOCC13),  I Y < 1 8) , ITM 02) , ISTAC20,  12), IVC4),  DELTAC  15) 

CQMMON/POTV'/P  C 20) , DEL, A 

COMMON/DIM/B, IBR 

COMMON/P 1 /NA, JK,  KZ 

COMNON/GRD/NPSIG,  NTF 

DATA  TSTCl), TST<2)/3HTST,  4H  SRC/ 


JK=4 
NA*  1 

CALL  STIND< IE,  2237,  0) 

URITEC4, 601) 

READC4,  6O0) A 
READC4,  600 ) B 
IBR=IFIXC10QGQ. *B) 

VJRITEC4,  2051) 

FORMAT C 1 X, 25HSPEC IFY  NO  OF  LINES  LT  16) 

READC4, 6004 ) NL I NES 
FORMAT (12) 

DELTX=. 5/ C FLOAT  <NLINES>) 

DELTAC 1 >=DELTX 

MTF-NLINES-1 

DO  2050  NT=2, NTF 

DELTACNT)=DELTACNT-1)+DELTX 

CONTINUE 

DELTAC 1 ) = 1 . / 1 8 . 

DELTA < 2 > -2 . /18. 

DELTAC 14) =8 . 6/18. 

DELTAC3) =4 . /IS. 

DELTAC4>=5. /18. 

DELTAC 15)=8. 7/18. 

DELTAC5 ) =7 . /18. 

DELTAC  6 ) = 7 . 3/18. 

DELTAC? ) =7 . 6/18. 

DELTAC  8 ) =3 . /IS. 

DELTAC9)=8. 1/18. 

DELTAC 10) =8. 2/18. 

DELTAC 1 1 )=8. 3/ 18. 

DELTAC 12>=3. 4/18. 

DELTAC 1 3>=8 . 5/13. 

FORMAT C 1 X/  ' INPUT  A, B, DELI,  DEL2, DEL3,  DEL4,  DEL5  : F5. 2,  5F5.  4' ) 
CONTINUE 

FORMAT C 1 X, 23HINPUT-  A, DEL:  F5. 2, F5.  4) 

URITEC4,  1 1 > DELTA <NA> 

DEL=DELTACNA) 

FORMAT < IX,  'DEL*',  F10.  4) 

CR=2. *BEL 
C 4 = CR*  r H3,  /6.  > 

DX1 l=A/6. 

BX21=A/6. 

DX3 1 = ( A-r C 6.  *C 4^/12. 

DX41=.CC4/2.  ) + U3>*A/6.  )-C4)/2. 

DX1 2=h/6. 

DX22=h/6. 

BX32=h/6. 
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DX42=C4 

DX52S  < h/2 . >-C4 
Cl=BX72''BX-2 
c2=ro:32  dxse 

602  FORMAT  < 1>£,  6<1X, F10.  4>) 

p<  n = i. 

P<2)-1 . 

P<3>-  01/<BN11*DX12> 

P < 4 >=S I NC3.  14159/6.  > 

P<5)=. 01/<DX11*DX22> 

P<6>- 1 . 

PC7)=SIN<3. 14159+2. /6. ) 

P < 8>  c 1 . 

P<9>=. 01 /CDX22+DX2 1 ) 

P<10>=. 0 1 / CDX32+DX21 ) 

PC l 1 >=SIN<3. 14159+3. /6. ) 

P < 1 2>  = . 0 1 / C DX32+DX3 1 ) 

P < 1 3 ) = 1 . 

P<14)=DX32/DX42 

PC  15)  = . 01/<DX31*DX42*C1) 

PC16)=DX32/DX52 

P< 17>=SIN<3.  14159+C1 . +CR>/2.  > 

PC 1S>=. 0 1 / < DX4 1 +DX42+C 1 > 

P < 1 9 ) =DX32/DX42 
P <20  > = . 01/<DX41+DX52*C2) 

IF<P<14>. LT. 1. 5)G0  TO  693 
PT0T1=PC14>+P<15) 

PC  14>  = 1 . 

P < 15 )=PTOT 1 
PT0T2=PC19)*P<18) 

P < 1 9 > * 1 . 

PC 18>=PT0T2 

698  CONTINUE 

IFCPC16). LT. 1. 5)G0  TO  699 
PTOT3=PC16)*P<20) 

PC 16>*1 . 

PC2G)=PT0T3 

699  CONTINUE 

DO  7G0  NP=1 / 20 
C WRITE < 4,  60 1 0) PC NP) 

700  CONTINUE 

600  FORMAT CF5. 2,  F5. 4) 

6010  FORMAT  C IX/ F10. 4) 

CALL  PDE 
CALL  MCON 
CALL  PDE2 
NA=NA+ 1 
JK= JK+ 1 

IFCHA. GT. NTF)GO  TO  3000 
GO  TO  900 
3000  CONTINUE 

CALL  DISK 

2021  FORMAT < IX*  T4,  ' ITM',  T14,  'XI',  T22,  'X2',  T3Q,  'X3',  T33, 
1 * X4  * » T4 6,  'X5', T54,  'X6', T62,  'XT',  T70,  'X8'  > 

C SEARCH  FOR  SPECIFIED  PSI  FOP  EOUIPOT  PRINTOUT 

4003  CONTINUE 

WPITEC4, 2006  > 
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2006  FORMAT < IX,  11 HSPECIFY  PSI> 

READ <4,  1 009  >FS I 
REAP<4, 1 009>P$ I D 
REAIK4,  1 021  >NPS  I 
1021  FORMAT < 1 2) 

HPSIG*0 

4002  COHTIIIUE 

DO  4000  IZR=I,NPSI 
1 PS I B I F I X < PS 1 * 1 00 . ) 

2022  FORMAT < IX, T9,  'P$l'* T20,  ' Y ' , T31,  'X',  T42,  'XI') 

K=  1 

NTF3=NTF+3 
DO  1000  H= 1 * HTF3 
DO  1001  NA  = 1,KZ 

I F < I ST  A < N, N A > . LT . I PS I > GO  TO  1002 

1001  CONTINUE 

1002  NB=NA 
NC=NB- 1 

DELSTA=FLOAT<ISTA<N,  NC)-ISTA(N,  NB>> 

DELTM=FlvOAT<ITM<NC)-ITH<NB)) 

Y< K>  = < FLOAT ( I TM<NC))-<DELTN* < (FLOAT C I STAC N,  NO -I PS I )^DELSTA> ) ) )✓ 
$ 10000. 

X=XLQC(N) 

XI=A-XLOC<N) 

K=K+ 1 

1000  CONTINUE 

1009  FORMAT  <2F1 0. 3,  12) 

1010  FORMAT < IX, 4 < 1 X, F10. 3)) 

IFCIZR. GT. 1 >G0  TO  4001 
IF<NPSIG. GT. 0>GO  TO  4001 
CALL  D.RU 

4001  CONTINUE 

CALL  DRWA 
PSI=PSI + PSID 
4000  CONTINUE 

PSI=PSI- (FLOAT (NPSI  )*PSID> 

NPSI G= 1 
WRITEC4,  1011) 

URITEC4,  1012) 

1011  FORMAT ( IX,  1 8HN0  GRID  OR  RESTART) 

1012  FORMAT (IX, 18HTYPE  2 FOR  RESTART) 

READC4,  102DMST 

IF (MST. EQ. 2) GO  TO  4003 

GO  TO  4O02 

STOP 

END 
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SUBROUTINE  PTE 

C PROGRAM  PDE**122673** 

DIMENSION  IPTC20)/ IPTVC20)/ ADELC IS) 

COMMON/Pl/NA,  JK/  K2 
COMMON/ POTV/P <20 )/  DEL/  A 

C0MH0N/DRUP/XL0C<18>,  IY<18>/  ITMC12)/  ISTA<20/  12)/  IVC4)/ DELTA< 15> 
DATA  IPT/2224/  2223/  2232/ 2235/ 2245/  2236/  2239/  2226/  2227/  2233/ 

1 2242/  2243/  2244/  2247/  2253/  2256/  2264/  2265/  2266/  2276/ 

CALL  LEX'CIE/  1) 

CALL  TSCALCIE/ 0) 

CALL  LOAD< IE) 

600  FORMAT CF5. 2/  5F5  4) 

405  CONTINUE 

CALL  LEX< IE/ 1) 

DEL=DELTACNA) 

XLOC< 1 >=A/6. 

XL0C<2>=2. *A/6. 

XL0C<3)=3. *A/6. 

ADELC  JK)*=DEL 

XLOCC JK)=3. *A/6. + DEL 

DO  750  IPV=1/ 20 

IPTV<IPV>  = IFIX< 10000.  *PCIPV)) 

750  CONTINUE 

CALL  INITACIE/0) 

CALL  CONSOCIE/0) 

CALL  LEXCIE/ 1) 

CALL  TSCALCIE/ 0) 

CALL  LOADC I E> 

5 CONTINUE 

CALL  STINDCIE/ 2277/  10G00> 

CALL  STINDCIE/ 2275/  0) 

DO  10  K= 1 / 20 

CALL  STINDCIE/ IPTCK)/ IPTVCK)) 

10  CONTINUE 

C SET  TIME  BASE 

CALL  STINDCIE/  2273/  1000) 

CALL  READC IE/  0200/  IDUM> 

CALL  LOADC IE) 

C URITEC4/  2000) 

2000  FORMAT C IX/ 27HSET  IC  POTS  260/ 261 , 262/ 263) 

RETURN 

END 
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SUBROUTINE  MCON 
INTEGER  PS  I < 100 ) 

COMMON  IJ, IK/ II J/ IDELX 
CALL  INITACIE/O) 

CALL  CCNSOCIE, 0) 

CALL  TSCALCIE/2) 

K=1 

200  I J=2225 
3K=2234 
JL— 2246 
I M=22?4 
11 J=0201 

IFCK. GT  1 )G0  TO  206 
IX=5000 

CALL  ST3N3KIE/  IJ/  IX) 

CALL  ST I ND < IE/  IK/  IX) 

CALL  ST I NIX  IE/  IL/  IX) 

CALL  STINIXIE/  IM,  IX) 

206  CALL  CONCIX/ PSI/ 1/ J> 

LX1=IX 
3J=2234 
I I J=U22 1 

201  IF <K. ED,  i)GO  TO  207 

GO  TO  208 

2B7  1X«5000 

GD  TD  20S 
208  3X=LX2 

2BS  CALL  CONCIX/PSI/  1/ J) 

LX2=IX 
IJ=2246 
II J=6241 

202  IFCK.EQ. 1)G0  TO  210 

GO  TO  211 

210  IX*5006 
GO  TO  212 

211  1 X=LX3 

212  CALL  CONCIX,  PSI/  1/ J) 

LX3=IX 

JJ=2274 
1 1 J=026 1 

203  IFcK. EQ.  1 > G 0 TO  213 

GO  TO  214 

213  1 X=5GO0 
GO  TO  215 

214  JX*LX4 

215  CALL  CONCIX/ PSI/  1/ J) 

LX4=3X 

K'=K+1 
3 J=2225 
33 J=&201 
3 X- LX 1 

CALL  REhLSI  <IX/ PSI/  1/ J/  IB) 

1 F CPS  I C I > GE. -100. AND. PSI <I). LE. 1O0)GO  TO  220 
GO  TO -226 
220  1J=2234 

II J=0221 
I X-LX2 
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CALL  READSKIX,  PSI,  I,  J,  IB> 

1F<F*SI<I>.GE.-100.  AND.  PSI < I > . LE.  100)  GO  TO  225 
GO  TO  203 

225  I J = 2 2 4 6 
II J=0241 
IX=LX3 

CALL  READS  I < IX,  PSI,  I,  J,  I B> 

IF<PSI<I>.  GE. -100. AND. PSI < I ) . LE.  100) GO  TO  235 
GO  TO  211 

226  GO  TO  206 

235  CONTINUE 

C PAUSE 

CALL  IC< IE> 

CALL  STIND<IE, 2277, 0> 

CALL  READ< IE, 0200, IVDUM) 

CALL  READCIE, 2222, IVDUN) 

CALL  UAIT <200) 

RETURN 

C STOP 

END 
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c 

c 

c 

c 

c 

c 

c 

5 


50 


15 


20 

21 


51 


25 

30 


52 

31 


53 


100 
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SUBROUTINE  CON< IX/ PS  1/  1/  J) 

INTEGER  P3IC100) 

COMMON  I J,  IK/  IIJz  IDELX 
CALL  READC IE>  2225/  1X325) 

CALL  READCIE/ 2225/ 1X325) 

CALL  WA I T < 70) 

CALL  READCIE/ 2210/  1X310) 

CALL  READCIE/ 2210/  1X310) 

CALL  WAIT C70) 

CALL  WAIT (70) 

1 = 1 

CALL  READSKIX/ PSI/  1/ J/  IB) 

IF< I . LE  1 )GO  TO  50 

IFCPSKI).  EQ.  PSICJ))GO  TO  900 

IFCPSKI).  GE. -100.  AND.  PSI  < I ) . LE.  100)  GO  TO  999 

IF< I . GT. l)GO  TO  15 

IFCPSICI).GT. 100) GO  TO  20 

GO  TO  100 

IFCI. EQ. l)GO  TO  20 

IF <P$I < I ) . LT. 0)GO  TO  999 

IF<PSI<I>. LT. PSI < J) )GO  TO  20 

GO  TO  100 

I X=I X+ I DELX 

1 = 1 + 1 

J=I-1 

CALL  READSKIX/ PSI/  1/ J/  IB) 

IF< I . LE. 1 )GO  TO  51 
IF<PSI<I). EQ. PSI < J) )GO  TO  900 
IFCPSKI).  GE. -100.  AND.  PSI<  I > . LE.  100)  GO  TO  999 
IFCPSKI).  LE.  -100) GO  TO*  999 
IFCPSKI  ).  GT.  PSICJ))GO  TO  25 
GO  TO  20 
IX=IX-IDELX 
1 = 1 

CALL  READSKIX/ PSI/  1/ J/  IB) 

IFCI. LE.  1 )GO  TO  52 

IFCPSKI).  EQ.  PSICJ))GO  TO  900 

IFCPSKI).  GE. -100.  AND.  PSIC  I ) . LE.  100)  GO  TO  999 

IX=IX-IDELX 

1 = 1 + 1 

J=I-1 

CALL  READSKIX/ PSI/  1/ J/  IB) 

IFCI. LE, l)GO  TO  53 
IFCPSKI).  EQ.  PSICJ))GO  TO  900 
IFCPSKI).  GE.  -100.  AND  PSIC  I ) . LE.  100)GO  TO  999 
IFCPSKI). LE.-1O0)GO  TO  999 
IFCPSKI  ).  GT.  PSICJ))GO  TO  20 
GO  TO  31 
I X= I X- IDELX 
1 = 1 + 1 
J=I-1 

CALL  READSKIX/ PSI/  1/ J/  IB) 

IFCI  LE. 1 )GO  TO  54 

IFCPSKI).  EQ.  PSI(  J)  )GO  TO  9O0 

IFCPSKI).  GE. -100.  AND  PSIC  I ) . LE.  1G0)GO  TO  999 

IFCPSKI).  GE.  1O0)GO  TO  999 

IFCPSKI).  GT.  PSKJ))GO  TO  180 
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110 


55 

111 


56 


900 

901 
999 


IX=IX+II»ELX 
1 = 1 

CALL  READSKIX.  RSI.  1,  J,  IB> 

IFC 1 . LE. 1 >G0  TO  55 
IF(PSI(I).E0  PSI(J))GO  TO  9O0 

IF(PSI(I). GE.-100.  AND,  PSICI). LE.  100>GO  TO  999 
IX=IX+ IDELX 
1 = 1 + 1 
J=I-1 

CALL  READSICIX,  PSI,  I,  J,  IB> 

I F < I . LE.  1 )GO  TO  56 

IF<PSI < I > . EQ. PS I < J> >GO  TO  900 

IF<PSI<I>.  GE.-100.  AND.  PSKD.LE.  100)GO  TO  999 
IF<PSK I >.  GE.  100) GO  TO  999 
IF <PSI <3>.LT.PSI < J>>GO  TO  100 
GO  TO  111 
WRITE < 4,  901) 

FORMAT < IX,  2HFU) 

CONTINUE 

RETURN 

END 
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SUBROUTINE  READSIOX,  PSI,  I, J, IB) 
INTEGER  PSI < 1 00) 

COMMON  IJ,  IK,  IIJ,  IDELX 
COMMONS  I M/B,  IBR 
I B=9000 
CALL  IC< IE) 

1006  FORMAT  OX,  3110,  2X,  110, 2X, 110) 

CALL  UA I T < 70 ) 

CALL  STINDOE,  IJ,  IX) 

C 

CALL  STINDOE,  2237,  IBR) 

CALL  WAIT (10) 

CALL  READOE,  0200,  ID2) 

CALL  WAITC70) 

C 

500  CONTINUE 

ANALOG  CONTROL  LOOP 
USES  ANALOG  COMPARATOR,  331 
CALL  0 POE) 

CALL  WAITC 1000) 

115  CALL  HOLDOE) 

CALL  UAIT<70) 

CALL  READOE,  IIJ,  IPSI) 

CALL  WA I T < 70 ) 

C 

CALL  ICOE) 

CALL  WAIT<30) 

CALL  STINDOE,  2237,  0) 

PSKI)*IPSI 
CALL  UAIT<70) 

CALL  READOE,  I J,  IXP) 

CALL  WAIT (70) 

C WRITE <4,  10G6)I  J,  IX,  IXP‘,  IPSI,  PSKI) 

CALL  WAITO00) 

I DELX= 1 0 

IFOA6SCPSIO)).  GE.  4000) GO  TO  10 
IF O ABS<PSI O ) ) . GE. 2000) GO  TO  9 
IFOABS<PSIO)).GE.  1000)GO  TO  8 
IF < I ABSCPSI < I ) ) . GE. 580) GO  TO  7 
IFOABSCPSICI)).  GE.  350)G0  TO  6 
I DELX= I DELX 
GO  TO  125 

6 I DELX=2* IDELX 
GO  TO  125 

7 I DELX=3* I DELX 
GO  TO  125 

8 IDELX=6*IDELX 
GO  TO  125 

9 IDELX=2* I DELX 
GO  TO  125 

10  I DELX- 14*1 DELX 

125  RETURN 

END 
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SUBROUTINE  PDE2 
C PROGRAM  PDE**1226?3** 

DIMENSION  I PT  <20 ) , IPTVC20),  ADELC1S) 

C0MM0N/P1/NA,  JK,  K2 
COMMON/POTV/PC20), DEL, A 

C0MM0N/DRUP/XL0CC18),  IYC1S),  ITMC12),  ISTAC20,  12),  IVC4),  DELTAC15) 
COMMON/D I M/B,  IBR 

DATA  IPT/2224,  2223,  2232,  2235, 2245,  2236, 2230, 2226,  2227,  2233, 

1 2242, 2243, 2244,  2247,  2253,  2256,  2264,  2265,  2266,  2276/ 

600  FORMAT CF5, 2,  5F5.  4> 

405  CONTINUE 

750  CONTINUE 

CALL  INITA< IE,  0) 

CALL  CONSOCIE,  0) 

CALL  LEXC IE, 0) 

CALL  ICC  IE) 

5 CONTINUE 

CALL  STINDCIE,  2277,  10000) 

CALL  STINDCIE,  2275,  O) 

C SET  TIME  BASE 

CALL  STINDCIE,  2273,  1000) 

ITMC1)*0 
CALL  WAITC70) 

C CALL  STINDCIE, 2275,  10000) 

CALL  UAITC70) 

K = 1 
MR=  1 

300  CONTINUE 

DO  3O0G  K=1 , 12 
Y=B*FL0ATCK)/1 1 
IYAS«IFIXC 10000.  *Y) 

CALL  WAIT <70 ) 

CALL  HOLDCIE) 

CALL  STINDCIE,  2237,  IYAS) 

CALL  WAIT C70) 

CALL  READCIE, 0200,  IVDUM) 

CALL  WAITC10O) 

CALL  READCIE, 0241,  ISTAC3, K)) 

CALL  WAIT C70) 

CALL  READCIE, 0221,  ISTAC2, K)> 

CALL  WAIT <70> 

CALL  READCIE, 0201,  ISTAC1, K)) 

CALL  WAITC70) 

CALL  READCIE, 0261,  ISTACJK,  K)> 

CALL  WAITC70) 

CALL  READCIE, 0271,  ITMCK)) 

CALL  WAITC70) 

C IFCITMCK). GE  IBR)GO  TO  102 

J1 =ISTAC 1, K) 

J2= I STAC  2,  K) 

J3= I ST  A C 3,  K) 

J4  = I STA  C4, K) 

I X- ITMCK) 

CALL  WAITC100) 

400  CONTINUE 

IFCK. EQ . 1 2 )G0  TO  102 
CALL  OPCIE) 
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c 

c 

c 

c 

c 

3000 

102 


200 


229 

228 

227 

231 


2001 


2005 

500 

2020 

2021 


CALL  WA 1 T C 1 900 > 

CALL  H0LIKIE) 

GO  TO  400 
I F (MR.  EC!.  1 >K=0 
MR=MR+ 1 

IF(K. GE. 12)G0  TO  102 

K=K  + 1 

GO  TO  30O 

CONTINUE 

CONTINUE 

CALL  UA IT < 100 > 

KZ=K 

CALL  IC(IE) 

CALL  UAITC200) 

FORMAT  < IX, 5(lX, 17)) 

CALL  IC(I€) 

CALL  WAITC1OO0) 

CALL  STINDCIE, 2275, 0) 

CALL  UAITC100) 

1=2225 

DO  2001  NI  = 1,  4 

GO  TO  (231,  227,  228,  229),  NI 

1=2274 

GO  TO  231 

1=2246 

GO  TO  231 

1=2234 

CALL  UAITC 1 00) 

CALL  READ<IE, I, IV(NI>) 

CALL  RE AD (IE, I, IV<NI>> 

CALL  WAIT <70 ) 

CONTINUE 
CALL  WAIT (70 ) 

URITE<4, 20O5) I VC1),  IV<2),  IVC3), IV(4) 

FORMAT (IX, 21HP0TS  225,  234, 246, 274, , 4C1X, 17)) 

CALL  WAIT (70) 

CONTINUE 

FORMAT  < 1 X, 9 ( 1 X,  I7>) 

FORMAT  ( IX,  T4,  'ITM',T14,  'X1',T22,  'X2',T30,  'X3',T38, 

X4 ' , T 46,  * X5 * , T54,  'X6',  T62,  'X7',T70,  'X8D 

RETURN 

END 
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SUBROUTINE  DISK 
DIMENSION  TST<2) 

COMMOH/UJ/ Y CIS)/  I PS  I 

COMHON/DRUP/XLOCOB),  IYU8), ITM02),  ISTAC20,  12),  IV<4),  DELTA < 15) 

CONMON/POTV/P<20),  DEL,  A 

COMMON/DIM/B,  IBR 

COHNON/Pl/NA,  JK,  K2 

COMMON/GRD/NPSIG,  NTF 

DATA  TST<i), TST<2)/3HTST,  4H  SRC/ 

CALL  ENTERC7, TST) 

NTF3=NTF +3 
DO  500  M=l, NTF3 
DO  500  NZ  = 1, KZ 

URITE<7, 2020) ITM(NZ),  ISTACM, NZ) 

500  CONTINUE 

2020  FORMAT OX,  SOX,  17)) 

CALL  CL0SEC7) 

RETURN 

END 
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SUBROUTINE  DRU 
COMMON/'W/'Y < 1 8 ) , IPSI 
CQMH0N/DR/Z<3S), ZYC38) 

COnMON/DRUP/XLOC<;S),  IY<1S>,  ITM<12>,  ISTAC20, 12), I V<4>,  DELTAC 15 ) 

C0MM0N/P0TWP<20>,DEL,  A 

COMMON/DItVB,  IBR 

COMHON/GRD/NPSIG, NTF 

NTF3=NTF +3 

DO  99  N=1 , NTF 3 

C R£AD<4,  9S)XL0CCN),  Y<N> 

B9  CONTINUE 

98  FORMAT  <2F10. 5) 

CALL  INITT(O) 

CALL  ERASE 

CALL  MDVABSC100, 100) 

CALL  DRWABS <100, 700) 

CALL  DRtJABS<  1000,  700) 

CALL  BRUABSC 1 000, 100) 

CALL  DRWABSC 100, 100) 

NLY=48 
LY=1 00 

NLYT=3FIX<10. *B)+48 
DO  251  N=l, 10 
CALL  MOVABSO00,  LY) 

CALL  BRUABSO0,  LY) 

CALL  MOVABSC50, LY) 

CALL  ANCH0C48) 

CALL  ANCH0C46) 

CALL  ANCHDCNLY) 

NLY=NLY+1 

IF  <NLY. GT. HLYT>GO  TO  261 
LY=<£00/<NLYT-48) )+LY 

251  CONTINUE 

261  CONTINUE 

DO  200  MT=1 , NTF3 
X3=A— XLDC  <t1T  ) 

KL=lF3X<XLDC<nT)*<900. /A>>+100 
KL1=1F1X<XI*<900.  /A>)+100 
CALL  MOVABSCKL, 90) 

CALL  DRUhBSCKL, 700) 

CALL  T1DVABS<KLI,  700) 

CALL  DRWABS<KLI, 90) 

CALL  MOVABS< KL”1 0/  80) 

XLDCX*XL0C<MT) 

1D1-IF3XC  XLDCX*1 O . ) 

ID2=IF3X<XLGCX*1 00. )-<10*IDl> 

XL0CX=XL0CX+A/8. 

I XC2=48 

I XC 1=4S 

DO  252  NR= 1 , 9 

IFCID1. EQ. HP ) IXC1=IXC1+NR 

1FCID2. EG. NR) I XC2=I XC2+NR 

252  CONTINUE 

CALL  ANCM0C4S) 

CALL  ANCHO < 1 XC1 ) 

CALL  AN CKO < I XC2> 
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253  CONTINUE 

CALL  M0VABS<20>  700> 
CALL  ANCH0C66) 

CALL  MOVABSC 1080/  30) 
CALL  ANCH0<65> 

200  CONTINUE 

CALL  M0VABS<50,  50) 

CALL  HOME 

CALL  AHMODE 

RETURN 

END 
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SUBROUTINE  DRWA 
COMMGN/W/YOS),  IPSI 
CQMH0N/DR/Z<36>,  ZY<36> 

CGMMON/DRWP/XLQC < 18)  ? IY<18),  ITM<12),  ISTA<20, 12), I V<4>, DELTA< 15) 

CQMMON/POTV/P<20>, DEL,  A 

COMMON/DIN/B, IBR 

COMMON /GRD/NPS I G,  NTF 

IF <NPSIG. NE.  1 )GO  TO  1000 

CALL  MOVABSC1 00,  100) 

CALL  DRUABSO00,  700) 

CALL  DRWABSaO00,  700) 

CALL  DRUABSC1000,  100) 

CALL  DRWABSC100,  10O) 

1000  CONTINUE 

NTF3=NTF+3 
DO  100  N=l, NTF3 
Z<N)=XLOC<N) 

Z<N+NTF3)=tt-XL0C<N> 

ZY  <N>  = Y < N) 

ZY<N+NTF3)=SYCN) 

100  CONTINUE 

I TOT= 1 

300  CONTINUE 

IF<ITOT. GT. 2000)GO  TO  400 
NTFR=  <2*NTF3)-1 
DO  220  N- 1 , NTFR 
IF<Z<N+1). LT. Z<N) )GO  TO  598 
GO  TO  220 

598  ZV1=Z(N+1 > 

ZV2=Z  <N ) 

ZY1=ZY<N+1) 

ZY2=ZY(N) 

Z<N)=ZV 1 
Z<N+1 )=ZV2 
ZYCN)=2Y1 
ZY<N+ 1 )=ZY2 
N=1 

ITOT=ITOT+l 
GO  TO  30O 
220  CONTINUE 

400  CONTINUE 

301  FORMAT <1X,  T5,  'X', T15,  ' Y' ) 

302  FORMAT < IX, 2F10.  5) 

CALL  HOME 

PSI=FLOAT  < IPSI )/l 00. 

X=0. 

DO  498  NP=1, 10O0 

TPS  1 = 1 00  *SIN<<3.  14159*X)/A) 

IFCTPSI.  GE.PSDGO  TO  497 
X=X+. 003 
498  CONTINUE 

497  IZX=IFIX<X*<9O0. /A))+10O 

XEND=A-X 

CALL  MOVABSCIZX,  100) 

NTFRA=NTFR*>1 
DO  411  NQ= 1 , NTFR A 
IF<NQ.  EO.  1 ) GO  TO*  473 
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473 


C 

465 


413 


414 


411 


C 


ZDEL*Z<NQ>-Z<NQ-1> 

IF<ZDEL. LT.  C. 001 >>G0  TO  411 
CONTINUE 

KLX-IFIX<2<NQ>*<900.  /A>>+100 
IFCKLX. LT. IZX>G0  TO  411 

IZXE=IFIX<XEND*<900.  /A>>* 100 
IF <KLX. GT.  1 2XE)G0  TO  411 
KLY=IFIX<ZY<NQ>*<60O. /B))+108 
I FCNQ . EQ. 9)G0  TO  413 
CONTINUE 

CALL  DRUAB$(KLX*  KLY) 

GO  TO  411 

CONTINUE 

IB1=IPSI/1D00 

ID2=<IPSI-CID1*1000>>/100 

I CX2S48 

I CX1*48 

DO  414  N-l# 9 

IF  < ID1 . EQ. N)ICX1  = ICX1+N 

IF< ID2. EQ. N>ICX2=ICX2+N 

CONTINUE 

CALL  ANCHOC ICX1 ) 

CALL  ANCHO< ICX2> 

CALL  MOVREL<-20,  0> 

GO  TO  465 
CONTINUE 

IZXA=IFIX<XEND*<900.  /A>)  + 100 
CALL  DRUABSCIZXA,  100> 

CALL  ANMODE 
STOP 
RETURN 
END 
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C-3.  Exact  Solution  Listings.  The  following  listings  are  used  for  the  exact  solution, 
which  is  run  using  “E IDEA”  since  the  exact  solution  also  required  chaining. 


$ £f  IT>Ef\ 

SUBROUTINE  POT  PflTf 

THIS  WILL  ACT  AS  THE  MAIN  PROGRAM 
DIMENSION  TST <2> 

CQMM0N/W/YC18),  IPSI 

C0MMQN/DRWP/XLQCU8),  I Y<  18),  ITM<  12),  ISTAC20,  12),  I V<4) , DELTA < 15) 
CQMMON/POTV/P<20), DEL, A 
COMMON/D I M/B, IBR 
COMMQN/P1/NA, JK, KZ 
COMMON/ GRB/NPS IG, NTF 
DATA  TSTC1),  TST <2) /3HTST,  4H  SRC/ 

JK=4 
NA=  1 

URITE<4,  601  ) 

REAB<4,  600) A 
REAB<4, 60O)B 
IBR=IFIX< 10000. *B) 

URI TE<4, 2051) 

2051  FORMAT < IX, 25HSPECIFY  NO  OF  LINES  LT  16) 

READ<4, 6004)NLINES 
6004  FORMAT  < 12) 

DELTX=. 5/ < FLOAT  <NLI NES) ) 

DELTA < 1 )=DELTX 
NTF =NLINES- 1 
DO  2050  NT=2, NTF 
DELTA<NT)=DELTA(NT-1 )+DELTX 
2050  CONTINUE 

DELThCI )=1. /18. 

DELTA<2)=2. /18 
DEL! A< 14)=8. 6/18. 

DELTA<3)=4. /IS. 

DELTA<4)=5. /18. 

DELTAC 15)=8. ?/18. 

DELTA <5 > =7. /18. 

DELTA  < 6 ) -7 . 3/18. 

DELTA<7)=7 . 6/18. 

DELTA<8)=8. /18. 

DELTA<9 )=8 . 1/18. 

DELTAC 10)=S. 2/18. 

BELT  A C 1 1 ) = 8 . 3/18. 

DELT A < 1 2 ).=8.  4/18. 

BELT  A < 1 3 ) = 8 . 5/18. 

601  FORMATS  IX,  * INPUT  A,  B,  DELI,  DEL2,  DEL3,  DEL4,  DELS  « F5. 2, 5F5. 4' ) 

900  CONTINUE 

697  FORMAT  C 1 X, 23H INPUT : A, DEL;  F5.2, F5.4) 

URITEC4, 1 1 )DELTA<NA) 

DEL=DELTACNA) 

11  FORMATOX,  * DEL=/ , F10.  4) 

CR=2 . *DEL 
C4=CR*A*<3.  /6 „ ) 

DX1 l=A/6. 

DX21=A/6. 

DX3 1 = ( A+ ( 6 . *C4))/12. 
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602 


698 


699 
C 

700 
60O 
6016 


3080 


DX4 1 =CC4/2.  ) + C C3.  *A/6.  )-C4>/2. 
DX12=A/6. 

DX22*A>6. 

DX32=A/6. 

DX42=C4 

DX52=CA/2. )-C4 

C1=I'X32/I.X42 

C2=DX32/DX52 

FORMAT  < IX/ 6C1X,  F10.  4>> 

PC1)=1. 

P C 2 ) = 1 

P<3>=. 01/CDX11*DX12> 

P<4)=SIN<3. 14159/6, > 

PC5)  = 01/<DXU*DX22> 

P<6>=1 . 

PC7)=SINC3. 14159*2. /6. > 

PC8>=1 . 

PC9)=. 01/C  DX22*DX2 1 ) 

P( 10)c. 01/CDX32*DX21 ) 

P< 1 1 >=SINC3. 14159*3. /6. ) 

PC  12)  = . 0I/<DX32*DX31> 

PC 13>=1 . 

PC14>=DX32/DX42 

P<  15)*=.  G1/CDX31*DX42*C1) 

PC16)=DX32/DX52 

PC 17>=SINC3.  141 59*  C 1 . +CR)/2.  > 

PC  1 8)tt. 01/CDX4 1*DX42*C1 ) 

PC 19>=DX32/DX42 
PC20>*  81/CDX4 1*DX52*C2) 
IFCPC14). LT.  1. 5)G0  TO  698 
PTOTl=PC 14>*PC 15) 

P C 1 4 ) « 1 

PC15)«PT0T1 

PT0T2=PC19)*PC13) 

PC 19)®  1 . 

PC 18>=PT0T2 
CONTINUE 

IFCPC16>. LT. 1 „ 5>GO  TO  699 
PTOT3=PC16)*PC20> 

PC 16)  = 1 . 

PC20>=PTOT3 
CONTINUE 
DO  700  NP= 1 > 20 
URITEC4, 6010  >PCNP> 

CONTINUE 

FORMAT  CF5. 2, F5. 4> 

FORMAT  C 1 X> F 1 0 . 4) 

CALL  EXACT 
NA=NA+ 1 
JK-JK+1 

IFCNA. GT. NTF)GO  TO  300O 

GO  TO  900 

CONTINUE 

NTF3eNTF+3 

DO  5O0  h* 1 t NTF3 

DO  500  NZ=1,KZ 

URITEC7>  2020) I TMCNZ)>  ISTACM, NZ) 
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500  CONTINUE 

202O  FORMAT  < IX/ 9<1X/  17)) 

2021  FORMAT  < IX/  T4,  ' ITM'/ T14/  'XI'/  T22/  'X2'/  T30/  'X3'/ T38/ 

1 'X4',  T46/  ' X5'  / T54/  'XS',  T62,  ' X7'  / T70/  'XS'  ) 

C SEARCH  FOR  SPECIFIED  PSI  FOR  EQUIPOT  PRINTOUT 

4003  CONTINUE 

WRITE  <4/ 2006) 

2006  FORMATOX/  1 1HSPECIFY  PSI) 

READC4, 1009)PS I 
READC4/ 1 009>PSID 
READ<4/  102DNPS1 
1021  FORMAT  < 12) 

NPS I G=0 

4002  CONTINUE 

DO  4000  I2R* 1 / NPSI 
I PS I = I F I X<  PS I ♦ 1 00.  ) 

2022  FORMAT < IX/ T9/  'PSI',  T20/  'Y'/ T31/  'X'/ T42/  .'XI'  > 

K=1 

NTF3=NTF +3 
DO  1000  N= 1 / NTF3 
DO  1001  NA=1/K2 

IF<ISTA<N/NA).LT  I PS I) GO  TO  1002 

1001  CONTINUE 

1002  NB=NA 
NC=NB-1  . 

DELSTA=FLOAT < ISTACN/  NO-ISTACN,  NB)  ) 

DELTM=FLOAT  <ITM<NC)-ITM<NB)) 

Y<K)=<FLOAT< ITM<NC))-<DELTM*<<FLOAT<ISTA<N/  NO-IPSI  )/DELSTA ) > ) ) / 
$ 10000, 

X=XL0C  < N ) 

XI*=A“XL0C<N) 

K=K+  1 

1000  CONTINUE 

1009  FORMAT <2F1 0.  3/  12) 

1010  F ORMAT < IX/  4<1X/F10.3>> 

IF < I2R. GT.  1 )G0  TO  4001 
IF<NPSIG. GT. 0)GO  TO  4001 
CALL  DRW 

4001  CONTINUE 

CALL  DRWA 
PSI=PSI+PSID 
4000  CONTINUE 

PSI =PSI-< FLOAT  <NPSI)*PSID) 

NPSIG=1 
WRITE<4/  101 1 ) 

WRITE  <4/  1012) 

1011  FORMAT < IX/  1 8HN0  GRID  OR  RESTART) 

1012  F0RMATC1X/ 1SHTYPE  2 FOR  RESTART) 

READC  4/  1021 )MST 

IF <MST . EG. 2) GO  TO  4003 

GO  TO  4002 

STOP 

END 
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SUBROUTINE  EXACT 

COMMON/DRWP/XLOC < 18)/  IY<13>,  ITM<12>,  ISTAC20,  12>,  IV<4>, 
1 DELTA< 1 5> 

COMMOH/POTV/P<20>, DEL, A 
COMMON/DIM/B, IBR 
COMMQM/P1/NA, JK, K2 
PI=3. 1415926 
XL0C<l>=A/6. 

XL0C<2>=2. *A/6. 

XL0C<3>=3.  *A ye. 

DO  5 1 = 4,  18 

XLOCC I >=XLOC < 3>+DELTA< I— 3> 

5 CONTINUE 

KZ=12 

DO  10  IX=1,  18 
DO  20  I YA=1, 12 
Y=B*FL0AT<IYA-1>/'11. 

ITM<IYA>=IFIX<Y*10000.  > 

Q1=PI*XL0C< IX) /A 

Q2=PI*B/A 

Q3=PI*<B-Y>/A 

PS 1 = 1 00.  *SIN<Q1>*<EXP<Q3>-EXP<-Q3>>/<EXP<Q2>-EXP<-Q2:>> 
I STA< I X,  IYA>  = IFIX< PS  1*100.  ) 

20  CONTINUE 

10  CONTINUE 

RETURN 
END 
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c 

99 

98 


251 

261 


252 


<nr\ 


SUBROUTINE  DRW 
COMMON/W/YCIS),  IPSI 
COMM0N/DR/ZC38), ZY<38> 

C0MM0!VDRUP/XL0C<18>,  IYC18),  ITM<12>,  ISTA<20,  12),  IV<4>,  DELTAC  15) 

COMMON/POTV/P<20>, DEL, A 

COMMON/DIM/B, IBR 

COMMOH/GRD/NPSIG, NTF 

NTF3=NTF+3 

DO  99  N=1 , NTF3 

READ<4,  9g>XL0C<N),  Y<N) 

CONTINUE 
FORMAT  C2F1G. 5) 

CALL  INITT <0) 

CALL  ERASE 

CALL  MOVABSU00,  100) 

CALL  DRWABS<100, 70O) 

CALL  DRUABSO000,  70O) 

CALL  DRWABS< 1 000, 100) 

CALL  DRUABS < 1 00,  100) 

NLY~4S 

LY=100 

NLYT^IFIXUO.  *B>+4 8 
DO  251  N= 1 , 10 
CALL  MOVABS< 100,  LY) 

CALL  DRWABSC90, LY) 

CALL  MOVABSC50, LY) 

CALL  ANCH0<43) 

CALL  ANCH0C46) 

CALL  ANCHO(NLY) 

NLY=NLY+1 

IFCNLY. GT. NLYT)GO  TO  261 
LY=<600/<NLYT-48) )+LY 
CONTINUE 
CONTINUE 

DO  2O0  MT-l, NTF3 
XI=A-XLOC<MT> 

KL=IFIX<XLOC<MT)*<900.  /A) >+100 
KLI=IFIX<XI *<900. /A) >+100 
CALL  MOVABSCKL, 90) 

CALL  DRWhBS<KL,  700) 

CALL  MOVABS<KLI, 70O) 

CALL  DRWABSCKLI, 90) 

CALL  MOVABSCKL-IO,  80) 

XLOCX=XLOC  <NT ) 

IDl=IFIX<XLOCX*l 0. ) 

ID2=IFIX<XL0CX+ 100. >-<10*IDl> 

XL0CX=XL0CX  + A/S . 

IXC2=48 

I XC1 =48 

DO  252  NR=1, 9 

I F< I D 1 . EQ. NR) IXC1=IXC1+NR 

IF< ID2. EQ. NR) IXC2- IXC2+NR 

CONTINUE 

CALL  ANCHORS) 

CALL  ANCHOdXCl) 

CALL  ANCHO < I XC2 ) 
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253 


200 


CONTINUE 

CALL  M0VABS<20,  700) 
CALL  ANCH0<*6> 

CALL  M0VABS<1000,  30) 
CALL  ANCHCK65) 
CONTINUE 

CALL  MOVABSC50,  50) 

CALL  HOME 

CALL  ANflODE 

RETURN 

END 
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•m  SC 


SUBROUTINE  DRW A 
COHMON/U/Y  < 1 8 > / I PS  I 
COMMON/ DR/ Z < 3S)/  ZYC35) 

C0NMGN/DRWP/XL0CC13),  I Y< IS) / ITMC12)/  I STAC 20/  12)/  IVC4)/  DELTAC 15> 

COMMON/ POTV/PC 20)/ DEL/ A 

COMMON/D  I M/B/  IBR 

COMMON/GRD/NPSIG/ NTF 

IF CNPSIG. NE.  1 )G0  TO  1OO0 

CALL  MOVABSCIQQ, 1Q0) 

CALL  DRUABSC100/ 700) 

CALL  DRWABSC1000/  700) 

CALL  DRUABSC1O00/ 100) 

CALL  DRWABSCIOO/  100) 

1000  CONTINUE 

NTF3=NTF+3 
DO  lOO  N* 1 / NTF3 
Z<N)=XLOCCN) 

ZCN+NTF3) *A-XLQCCN> 

ZY<N)=YCH) 

ZYCN+NTF3)=YCN) 

100  CONTINUE 

I TOT= 1 

300  CONTINUE 

IF< ITOT. GT. 2000) GO  TO  400 
NTFR=<2*NTF3)-1 
DO  220  N=l/ NTFR 
IFCZCN+l). LT. ZCN))GO  TO  59S 
GO  TO  220 

598  ZV1=Z<N+1 ) 

ZV2=Z(N) 

ZYl =2YCN+ 1 ) 

ZY2=ZYCN) 

ZCN) =ZV 1 
ZCN+l )=2V2 
2YCN)=ZY1 
ZYCN+ 1 )=ZY2 
N=  1 

IT0T=IT0T+1 
GO  TO  360 
220  CONTINUE 

400  CONTINUE 

301  FORMAT C 1 X/ T5/  'X'/ T15/  * Y' ) 

302  FORMATC IX/ 2F10. 5) 

CALL  HOME 

PSI=FLOAT< IPSI )/100. 

X=0. 

DO  498  NP- 1/ 1080 

TPSI=1O0. *$IN<<3. 14159*X)/A> 

IF < TPS I . GE. PS I) GO  TO  497 
X=X+ . 005 
498  CONTINUE 

497  I ZX- I F I X < X*  C 980 . /A))+100 

xeni-h-x 

CALL  MOVABSC IZX/  100) 

NTFRA=NTFR+1 
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473 

C 

465 

413 

414 

411 

C 


HO  411  NQ=1,NTFRA 
IFCNG. EC. 1 >G0  TO  473 
ZDEL*ZCHQ>-Z<N0-1> 

IFCZDEL.LT.  <.001>>G0  TO  411 

CONTINUE 

KLX* IF1XC2C HQ )*< 900 . ✓A) > + 100 
IFCKLX.LY  IZX)GO  TO  411 

IZXE*IFIXCXEND*<900.  <^>>+100 
I F CKLX* GT.  IZXE>CO  TO  411 
KLY=IFIXCZYCNQ)*<600.  /B>)  + 10O 
IFCHQ. EQ. 5>G0  TO  413 
CONTINUE 

CALL  DRWABSCKLX,  KLY> 

GO  TO  411 

CONTINUE 

ID1*IPSI/10O0 

ID2=<IPSI-aDl*lO00)  >/100 

I 0X2=48 

ICX1 *48 

DO  414  N»l,  9 

IFCIDl . EQ. N> ICX1=ICX1+N 

IFCID2. EQ. N> ICX2*ICX2+N 

CONTINUE 

CALL  ANCHOCICXl) 

CALL  ANCHOC ICX2) 

CALL  MOVRELC-20,  0> 

GO  TO  465 
CONTINUE 

IZXA*IFIXCXEND*<900.  sft))+100 
CALL  DRUABSCIZXA,  100> 

CALL  ANMODE 
STOP 
RETURN 
END 
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C-4.  Stored  Data  for  \jj(x,  y).  The  ^(y)-data  taken  for  each  x-station  during  the 
exact  and  hybrid  solutions  are  provided  as  comparison  data  between  solutions. 


3 

4993 

126 

4652 

364 

3990 

545 

3504 

728 

3032 

911 

2567 

1091 

2112 

1273 

1665 

1455 

122$ 

1636 

787 

1819 

356 

2000 

-71 

3 

8650 

126 

8048 

364 

6894 

545 

6043 

728 

5210 

911 

4390 

1091 

3584 

1273 

2788 

1455 

1999 

1636 

1208 

1819 

418 

2000 

-373 

3 

9987 

126 

9312 

364 

8014 

545 

7855 

728 

6116 

911 

5196 

1091 

4293 

1273 

3405 

1455 

2526 

1636 

1652 

1819 

788 

2000 

-72 

3 

9838 

126 

9192 

364 

7849 

545 

6878 

728 

5924 

911 

4985 

1091 

4057 

1273 

3137 

1455 

2222 

1 636 

1305 

1819 

383 

2O00 

-547 

3 

9386 

126 

8781 

364 

752S 

545 

6624 

728 

5736 

62 


i 


911 

4863 

1091 

4002 

1273 

3151 

1455 

2306 

1 636 

1466 

1819 

628 

20C0 

-21 1 

3 

7648 

126 

7072 

364 

6118 

545 

5372 

728 

4642 

911 

3925 

1091 

3222 

1273 

2528 

1455 

1838 

1636 

1153 

1819 

473 

2O00 

-209 

3 

6415 

126 

5940 

364 

5122 

545 

4495 

728 

3882 

911 

3280 

1091 

2690 

1273 

2110 

1455 

1533 

1636 

962 

1819 

395 

2000 

-175 

3 

3414 

126 

3160 

364 

2684 

545 

2338 

728 

2004 

911 

1679 

1091 

1 364 

1273 

1055 

1455 

750 

1636 

448 

1819 

143 

2000 

-156 

3 

2920 

126 

2734 

364 

2282 

545 

1984 

728 

1695 

911 

1415 

1091 

1144 

1273 

878 

1455 

614 

1636 

354 

1819 

96 

2O00 

-162 

3 

2415 

126 

2223 

364 

1874 

545 

1621 

728 

1 o 3 0 

91  1 

1143 

1091 

922 

1273 

702 

1455 

464 

1 636 

269 

1819 

52 

2000 

-165 

3 

1731 

126 

1588 

364 

1326 

545 

1140 

728 

963 

911 

796 

1091 

634 

1273 

475 

1455 

320 

1636 

164 

1819 

9 

2000 

-149 

3 

1562 

126 

1430 

364 

1 189 

545 

1022 

728 

860 

911 

705 

1091 

559 

1273 

416 

1455 

271 

1636 

134 

1819 

-7 

2000 

-152 

64 


(T  y / r T 


e 

5000 

181 

4495 

363 

40O4 

345 

3527 

727 

3061 

969 

2606 

1090 

2158 

1272 

1718 

1454 

1284 

1636 

853 

1818 

426 

2000 

0 

0 

8660 

181 

7785 

363 

6936 

545 

6109 

727 

5303 

909 

4513 

1090 

3739 

1272 

2976 

1454 

2224 

1636 

1478 

1818 

738 

2000 

0 

0 

10000 

181 

8990 

363 

8009 

545 

7055 

727 

6123 

909 

5212 

1090 

4317 

1272 

3437 

1454 

2568 

1636 

1707 

1818 

852 

2000 

0 

0 

9848 

181 

8853 

363 

7887 

545 

6947 

72? 

6038 

909 

5132 

1090 

4252 

1272 

3385 

1454 

2529 

1636 

1681 

1818 

839 

2000 

0 

0 

9396 

181 

8447 

363 

7526 

545 

6629 

72? 

5754 

909 

489? 

1G90 

405? 

1272 

3230 

65 


1454 

2413 

1636 

1604 

1818 

SO0 

2000 

0 

0 

7660 

1S1 

6886 

363 

6135 

545 

5404 

727 

4698 

909 

3992 

1090 

3307 

1272 

2633 

1454 

1967 

1636 

1308 

1818 

652 

2000 

0 

0 

6427 

181 

5778 

363 

5148 

545 

4534 

727 

3936 

909 

3358 

1090 

2775 

1272 

2209 

1454 

1650 

1636 

1097 

1818 

547 

2000 

0 

0 

3420 

181 

3074 

363 

2739 

545 

2412 

727 

2094 

909 

1782 

1690 

1476 

1272 

1175 

1454 

878 

1636 

534 

1818 

291 

2000 

0 

0 

2923 

181 
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545 
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66 


1090 

1044 

1272 
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1454 

621 

1636 

413 

1818 

206 

2000 

O 

0 

1736 

C-5.  Program  Control  Flow  Charts.  The  hybrid  program  is  shown  in  Figure  C-l, 
while  the  exact  program  is  depicted  in  Figure  C-2. 
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Cs 

CO 


SUB  PDE2  f- 


SUB  DRW 


SUB  DRWA 


MAIN 


CALL  PDE 


CALL  MCON 


CALL  PDE 2 


CALL  DISK 


CALL  DRW 


CALL  DRWA 


SUB  PDE 


SUB  CON 


CALL  READS  I 


SUB  MCON 


CALL  CON 


CALL  READS  I 


^ 


SUB  READS  I 


Figure  C-l.  Block  diagram  — hybrid  program. 


SUB  DISK  > 


SUB  DRW 


MAIN 
CALL  PDE 
CALL  EXACT 
CALL  DISK 
CALL  DRW 
CALL  DRWA 


Figure  C-2.  Block  diagram  — exact  program. 


SUB  PDE  | 


~SUB  EXACT  f 


^ SUB  DRWA 


C-6.  Chaining  Routine.  The  chaining  routine  is  as  follows: 
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define  reside:::  code 
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VA  I T 

3/356-34/06 
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rrcr 
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